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We analyze the collective behavior of a lattice model of 
pulse-coupled oscillators. By means of computer simulations 
we find the relation between the intrinsic dynamics of each 
member of the population and their mutual interaction that 
ensures, in a general context, the existence of a fully synchro- 
nized regime. This condition turns out to be the same than 
the obtained for the globally coupled population. When the 
condition is not completely satisfied we find different spatial 
structures. This also gives some hints about self-organized 
criticality. 



The collective behavior of large assemblies of pulse- 
coupled oscillators has been investigated quite often in 
the last years. Many physical and biological systems can 
be described in terms of populations of units that evolve 
in time according to a certain intrinsic dynamics and in- 
teract when they reach a threshold value ^ . Although it 
was known long ago that the members of these systems 
tend to have a synchronous temporal activity, a rigor- 
ous treatment of the problem has been considered only 
in the last decade |^|-^. Up to now, the most important 
efforts have been focused on systems with long-range in- 
teractions because in this case analytical results can be 
derived by applying a mean-field formalism. Relevant is 
the work by MiroUo and Strogatz (MS) Q who discovered 
under which conditions mutual synchronization emerges 
as the stationary configuration of the population. Later, 
the study has been generalized to other situations [p|"p^ . 

When the oscillators form a finite dimensional lattice 
where only short-range interactions are allowed, the spec- 
trum of behaviors is more complex. For instance, some 
recent papers have related lattice models of pulse-coupled 
oscillators to models displaying self-organized criticality 
(SOC) pll^, systems which self-organize, due to its 
own dynamics, into a critical state with no characteris- 
tic time or length scales Of particular interest for 
us is the model proposed by Feder and Feder (FF) ||l^] 
to study stick-slip processes in earthquake dynamics. In 
such model each cell of a 2d square lattice is described 
physically in terms of a state variable E{t) (hereafter 
called energy) that evolves linearly with time. Once the 
energy of a cell reaches a threshold value {Eij > Ec — 1) 
it becomes critical and "fires" transferring energy to its 
nearest neighboring cells according to the following rules 
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where e is the strength of the coupling. In turn, some of 
these neighboring cells may become critical generating an 
avalanche that propagates through the lattice. When one 
avalanche is triggered the intrinsic dynamics is stopped 
and only when it is over does the driving act again. In 
this way, there is a clear separation of two time scales. 
By identifying the state variable E with a voltagelike 
magnitude one can establish an analogy between the FF 
model and some models of integrate- and-fire oscillators. 

Without any other ingredient the FF model displays, 
in the stationary state, and for open boundary condi- 
tions, relaxation oscillations (RO) which give rise to spa- 
tial structures formed by large assemblies of units, all 
with the same phase. In this sense we could talk about a 
macroscopic (local) degree of synchronization. However, 
when a dynamical noise is added to the system new col- 
lective behaviors appear, since the synchronized state is 
no longer an attractor of the dynamics. For instance, for 
£ — 0.25 1 11 1, the distribution of avalanche sizes follows 
a power-law decay characteristic of SOC. 

In the short range models described above the en- 
ergy of each oscillator varies with a constant driving rate 
f{E) = ^ — C . The analysis has been extended to lin- 
ear driving rates finding new conditions to observe 
RO even in presence of noise. It would be interesting to 
consider a more general context, where f'{E) can take 
an arbitrary shape with the only restriction, f{E) > 0. 
There is another underlying hypothesis in (|l|) that re- 
stricts the range of interest of these models: the strength 
of the coupling e is constant. In many studies of biolog- 
ical pacemakers it is assumed that the response of a cell 
due to an external stimulus is a function of the phase, 
A(0), which depends on the current state of the oscilla- 
tors |l8| , |l9| . The phase-shift induces an energy-shift e{E) 
(hereafter called the energy response curve (ERG)) that, 
in general, is a non-constant function of E. Our goal 
is to investigate a wide variety of models with arbitrary 
ERG and driving rate f{E) beyond those considered in 
previous works. Although in such a general situation a 
broad range of different regimes can be observed, we will 
focus our study on the conditions required for the system 
to develop a fully synchronous stationary state, empha- 
sizing the relevance of the response of a given oscillator 
at its reset point, e{E — 0). The analysis also provides 
information about the possible nature of SOG. Our in- 
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terest mainly concerns coupled map lattice models, but 
to start the discussion a mean-field model is introduced 
to clarify several dynamic aspects. 

Let us consider a system of globally coupled oscillators 
which interact through a given s{E). The dynamics of 
each unit is given by 



dt 



f{E,)+e{E,)5{t-t,) 



(2) 



plus the reset condition for Ei > 1. Here tj denotes the 
time at which unit or group j fires. In this description 
both time scales (driving and coupling) appear in the 
same equation. We also could consider another equiv- 
alent description where both time scales are separated 
such as in (1), but this fact implies to substitute e{E) by 
£{E) related through 



E+eiE) 



siE') 



S{t - tj)dt = 1. 



(3) 



To go from e{E) to e{E) is trivial, but the inverse implies 
to deal with an integral equation that, in general, can 
only be solved numerically. Now, it is simple to derive 
sufficient conditions to ensure perfect entrainment be- 
tween both oscillators. The method consists of applying 
the following transformation, used in different contexts 
by several authors SJSQ], 



y = £o 



dE 







(4) 



where Eq is defined to ensure that y{l) = 1. By substi- 
tuting into (p|) we get 



dyi 
dt 



eo—^ + eoS{t-tj). 
e{Ei) 



(5) 



Now, the evolution of the system is described in terms 
of a new variable y for which the coupling is constant. 
A case of particular interest is that of a zero advance 
at the reset point {e{E = 0) = 0). In such case, this 
transformation is well defined if the energy transferred 
Vy ^ is constant (eo) except for y = which is exactly 
zero. This condition plays the role of a refractory time, 
which provokes that the units that have fired have zero 
phase at the end of the interaction process. 

A perfect analogy can be established with the prob- 
lem studied in MS They showed that for a constant 
positive coupling, no matter how small, the synchronized 
state is an absorbing state of the dynamics if the driv- 
ing rate is positive and its first derivative negative. By 
applying these conditions to the new function 
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we find the following relation between the driving rate 
and e{E) that ensure the synchronization of the popula- 
tion: 
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This means that given the features of the driving rate it 
is always possible to find the oscillators firing in unison 
provided one chooses the suitable ERC. 

The interesting point is to know whether the mathe- 
matical result we have derived above can be extended to 
networks with short-range interactions. We have consid- 
ered FF coupling (|^) with ejE) and nonconstant driving 
rate with dynamical noise ||21[| . It is important to realize 
that in contrast with mean-field models, where synchro- 
nization emerges in a process where clusters of oscillators 
of increasing size merge with each other (absorption pro- 
cess) and never break up, in a coupled map lattice model 
big assemblies of oscillators with the same phase (which 
eventually may break up) are generated through large RO 
that sweep the whole lattice (avalanches of the size of the 
system). Then, for constant, instantaneous couplings and 
due to the different nature of both mechanisms, the con- 
ditions required to find synchronization/RO are different 
p2[ . However, the situation may change if one consid- 
ers a non-constant ERC. To analyze this new situation 
it is convenient to apply the same arguments discussed 
in but now for a state-dependent coupling. For a 
square lattice with open boundary conditions, it is not 
difficult to show that a necessary condition to observe 
RO of the size of the system is 



E{1 - 0(4e(O))) + 2e(l - (/'(4e(0))) > 1 



(8) 



where the energy and the phase 
the following expression 
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are related through 
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First of all, we observe the relevance of the response 
of a cell at the reset value. In general, a non-zero ad- 
vance at this point (e(0) 7^ 0) gives rise to a spatial dis- 
tribution of phases after an avalanche; it is only under 
these conditions that SOC can be obtained |ll|. On the 
other hand, it is clear that if e{0) — the inequality (H) 
is always satisfied, for any driving rate and ERC. This 
means that the appearance of RO involving all the sites 
implies a perfect synchrony between all the elements of 
the lattice, boundaries included. However, (^) does not 
ensure the existence of those RO. It is a necessary but 
not a sufficient condition for the system to achieve a com- 
plete synchronization. We will check by simulations that 
condition (j^) is, as in the long-range case, sufficient to 
produce synchrony. 

Before giving evidence of this fact, we want to men- 
tion some conclusions that can be extracted from this 
result. Although, in general, one cannot map straightfor- 
wardly results from mean-field theories to coupled map 
lattices models, for these particular systems and as long 
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as synchronization is concerned, conditions (|7|) which are 
strictly derived in a mean-field frame, can be applied to 
short range systems. This assumption is also in agree- 
ment with a conjecture proposed in [Q for the special 
case of fiE) < and e\E) = 0. 

Several models have been considered in our simulations 
finding a complete agreement with the theory. In par- 
ticular, we have performed simulations for the Peskin's 
model |2^] , where the driving rate is given by 

f{E)^j{K~E), (10) 

where 7 gives the slope of the driving rate and K = 
(1 — e^'') . Starting from a random distribution of 
phases and for 7 > {f'{E) < 0) we have observed 
that the population always synchronizes when e{E) is 
a positive function of E with e'{E) > 0, provided that 
e(0) = 0. However, for 7 = a monotonously increasing 
e{E) is needed. 

FIG. 1. Number of avalanches (filled) and time in period 
units (hollow) that a 32 x 32 lattice needs to get a complete 
synchronization as a function of 7 — 7' in a log — log scale. 
The symbols are averages over 25 realizations and the error 
bars (which for the avalanches are of the size of the symbol) 
correspond to the standard deviation. The straight line shows 
the (7 — 7')"^ behavior. 

To complete this idea we have studied the time re- 
quired for the population to synchronize, starting with 
random phases between and 1, for an ERC given by 

e{E)=eoj'{K' ^E), (11) 

where K' — (^1 — e^"* ^ , in a 32 x 32 lattice. We have 
plotted our results in Fig. ^ for open boundary condi- 
tions. There we can see how this time and the number of 
avalanches diverge as the difference between 7 and 7' ap- 
proaches zero. We have also checked a driving rate and 
an ERC both given by power laws and the results are 
very similar; the time required to synchronize the system 
grows as (a — a')~^, where a and o! are the exponents 
of the driving rate and of the ERC, respectively. Finally 
we have simulated a driving rate of the type (^) with 
a power-law ERC. For an exponent a' > 1 the system 
synchronizes for any value of 7; in particular the time 
needed to obtain the fully synchronized state grows ex- 
ponentially as e~^/'^ when a' = 1. On the other hand, 
for a' < 1 there is only a range of values of 7 for which 
we get synchronization. 

All these results arc in agreement with our prediction 
(^. At this point it is important to remark that when 
the driving and the ERC have the same functional de- 
pendence on E the inequality can be satisfied for all E 
and the divergence of the time appears when both curves 
/'// and e'/e approach each other. However when the 
dependence is different both curves can cross and this 



leads to a more complex behavior. Finally, the case of a 
power-law ERC with a' > 1 is very important because in 
this case the transformation (Q) cannot be performed and 
therefore no mapping with the MS result can be done; 
nevertheless (^ still represents the sufficient condition to 
get a fully synchronized regime. Lattices with periodic 
boundary conditions have also been checked and the con- 
clusions are the same than for open ones; furthermore, 
the time required to synchronize scales in the same way. 

We have also investigated the behavior of the model 
when one does not expect synchronization for a constant 
ERC, provided that £(£■ = 0) = 0. First of all, let us 
recall the behavior of two coupled oscillators The 
phase of the one oscillator when the other one arrives to 
the threshold transforms according to 

00 ^ l-0(£;(0o)+e) (12) 

provided that E{(j)Q) + e < 1, otherwise the oscillators 
synchronize. This transformation has always at least 
one fixed point, (j)Q, which leads to different behaviors 
depending on the slope of f{E). Thus, for f'{E) < 
(Vi?), is unique and unstable and the two oscil- 
lators will always synchronize for any positive value of e 
@. However, for f'{E) > (VE) the stabifity of the 
fixed point changes, and (jj^ becomes an attractor, which 
means that the oscillators can either be phase-locked, 
when e < 1 — i?(0o)j or synchronized, otherwise. For the 
Peskin's model ( p^ 4>q corresponds to 

= 1/2 - l/7sinh"i (£sinh(7/2)) . (13) 

This result enables us to perform a qualitative anal- 
ysis of the lattice model with nearest-neighbor interac- 
tions. Starting from a random distribution of phases, 
and only for periodic boundary conditions, the oscilla- 
tors will tend to be locally synchronized, or phase-locked 
depending on their phases and the parameters' values. 
This makes us to suspect the existence of well defined 
spatiotemporal structures of phase-locked oscillators in 
the stationary state, of the same kind as those described 
in p3|] , for which simple return maps can be written. 
The most simple of these configurations one can imag- 
ine is a chessboard lattice where "black" sites have the 
same value of the phase (0o) when all the "white" sites 
arrive to the threshold. Once the white ones have fired, 
the black ones are driven up to the threshold value and 
one gets the same structure as before with a phase 00 
that transforms according to (^2|) replacing e by 4£. This 
means that there exists a fixed point for this structure, 
whose value is the same than that obtained from ( [T^ ) 
replacing e by 4e, that is an attractor for the dynamics 
of the lattice. It is important to remark that this struc- 
ture is characterized by the fact that a given oscillator is 
phase-locked with its four neighbors. Other simple struc- 
tures arise when an oscillator can be synchronized with 
one neighbor and phase-locked with the remaining ones 
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(all these neighbors with the same phase), and so on, 
as far as the number of synchronized neighbors is con- 
stant through the lattice. What we have observed in our 
simulations is that these structures are attractors of the 
dynamics within different domains of the space of ini- 
tial conditions. The corresponding fixed points are given 
by ( p^ ) changing e by 4e, 3e, or 2e. Nevertheless, these 
structures are not the only attractors since more compli- 
cated patterns involving more than two different phases 
also can be built. Among of all these configurations, the 
"chessboard" one is the most relevant because it has the 
larger domain of attraction, although the relative weight 
depends on the parameters as well as on the system size, 
and it is the most stable one in the sense that small fluc- 
tuations break the other structures in its favor We 
want to point out that these phase-locked states are char- 
acteristic of lattice models with short-range interactions, 
since they have no counterpart with models of all-to-all 
pulse-coupled oscillators. 

Finally, we believe that the tendency of our model ei- 
ther to synchronize or to form phase-locked structures al- 
lows us to go one step further in the current understand- 
ing of SOC phenomena. Middleton and Tang ||l^ noticed 
that SOC appears, in a uniformly driven model and a 
slightly different coupling, as a consequence of marginally 
stable phase- locking between neighbors. We have shown 
how this marginal stability (which corresponds to the 
equality in the r.h.s. of (|^)) is broken in favor of a syn- 
chronization or a stable phase-locking depending on the 
driving rate and on the ERC. Thus, although further in- 
vestigation along this line would be required, one can 
conjecture that SOC is critical in the sense that balances 
the tendency into one or another direction. 

In this paper we have shown how is possible to reduce 
a general population of biological oscillators to a simple 
model with constant ERC, that allows analytical results 
for all-to-all coupling by means of a very easy transfor- 
mation. Surprisingly, short-range simulations verify the 
same conditions for the synchronization that the long- 
range version of the model. Moreover, we also find new 
states of the system with no analogous in the long-range 
case. These behaviors allow us to give some hints about 
the origin of SOC. 

The authors are indebted to L.F. Abbott, K. Chris- 
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